
#include <imcmc/imcmc.hpp>
#include <imcmc/parser++.hpp>
#include "CosmoTheory.hpp"
#include "RSD.hpp"

using namespace std;
using namespace imcmc;
using namespace imcmc::parser;

void Compute_RSD_Vals(  CosmoTheory* cosmo,
                        double* zeff,
                        double* rsd,
                        int&    size){
    for( int i=0; i<size; ++i ){
        double f, sigma8;
        f       = cosmo->engine->get_f(zeff[i]);
        sigma8  = cosmo->engine->get_sigma8(zeff[i]);
        rsd[i]  = f*sigma8;

        // cout << "## z = " << zeff[i]
        //      << ", f = " << f
        //      << ", sigma8 = " << sigma8 << endl;
    }

    // exit(0);
}

double Likelihood_RSD(  imcmc::imcmc_double&    param,
                        double&                 lndet,
                        double&                 chisq,
                        void*                   model,
                        void*                   data,
                        imcmc::istate&          state ){

    lndet = chisq = 0.0;

    CosmoTheory	*cosmo  = static_cast<CosmoTheory*>(model);
    Data_RSD    *rsd 	= static_cast<Data_RSD*>(data);

    double *rsd_prediction = new double[rsd->size];
    Compute_RSD_Vals(	cosmo,
                        rsd->zeff,
                        rsd_prediction,
                        rsd->size );
    //  debug
    for( int i=0; i<rsd->size; ++i ) {
        // cout << "z = " << rsd->zeff[i]
        //      << ", Theory value " << rsd_prediction[i]
        //      << ", Obs value " << rsd->val[i] << endl;
        chisq += pow((rsd_prediction[i]-rsd->val[i])/rsd->err[i],2);
    }
    // exit(0);    

    delete[] rsd_prediction;

    return lndet - 0.5*chisq;
}